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Abstract 

Background: The serious feeding- and microbiota-associated intestinal disease, necrotizing enterocolitis (NEC), 
occurs mainly in infants born prematurely (5-10% of all newborns) and most frequently after formula-feeding. 
We hypothesized that changes in gene methylation is involved in the prenatal maturation of the intestine and 
its response to the first days of formula feeding, potentially leading to NEC in preterm pigs used as models for 
preterm infants. 

Results: Reduced Representation Bisulfite Sequencing (RRBS) was used to assess if changes in intestinal DNA 
methylation are associated with formula-induced NEC outbreak and advancing age from 10 days before birth to 
4 days after birth. Selected key genes with differentially methylated gene regions (DMRs) between groups were 
further validated by HiSeq-based bisulfite sequencing PCR and RT-qPCR to assess methylation and expression levels. 
Consistent with the maturation of many intestinal functions in the perinatal period, methylation level of most genes 
decreased with advancing pre- and postnatal age. The highest number of DMRs was identified between the newborn 
and 4 d-old preterm pigs. There were few intestinal DMR differences between unaffected pigs and pigs with initial 
evidence of NEC. In the 4 d-old formula-fed preterm pigs, four genes associated with intestinal metabolism (CYP2W1, 
GPR146, TOPI MT, CEND1) showed significant hyper-methylation in their promoter CGIs, and thus, down-regulated 
transcription. Methylation-driven down-regulation of such genes may predispose the immature intestine to later 
metabolic dysfunctions and severe NEC lesions. 

Conclusions: Pre- and postnatal changes in intestinal DNA methylation may contribute to high NEC sensitivity in 
preterm neonates. Optimizing gene methylation changes via environmental stimuli (e.g. diet, nutrition, gut microbiota), 
may help to make immature newborn infants more resistant to gut dysfunctions, both short and long term. 
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Background 

Epigenetics represent stable and heritable changes in 
gene expression without changing the DNA sequence, 
providing a mechanism whereby environmental factors 
can affect functions of specific cells, tissues and organs 
[1]. DNA methylation is the best known mechanism of 
epigenetic modification [2] and it is particularly dynamic 
during the embryonic phase of mammalian development 
[3]. A second phase of intense DNA methylation alterations 
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may take place during the perinatal period, when many 
organs and tissues undergo dramatic changes in both 
physiological and environmental conditions. Genome-wide 
epigenetic changes may also be induced by specific nutri- 
ents [4] and such methylation changes can temporarily or 
permanently alter the expression of tissue-specific genes, 
thereby supporting the short and long term maturation of 
tissues and organs. 

Around birth, the mammalian gastrointestinal tract 
(GIT) has to adapt rapidly to the change in nutrition 
mode (from placental to enteral nutrient intake) and 
living environment (from a sterile life in utero to extra- 
uterine exposure to bacteria and viruses) [5]. The exact 
time frame of this critical adaptation period varies widely 
among the different mammalian species [5] but in pigs 
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intestinal structure and function change markedly from 
ten days before birth to 3-4 days after birth. In pigs, 
preterm delivery at 88-95% gestation is associated with 
organ immaturities and clinical complications similar to 
those in infants born at 70-90% gestation [5]. Conse- 
quently, preterm infants and pigs lack the final matur- 
ation of the GIT and suffer from a series of intestinal 
defects, leading to increased risk of intestinal compli- 
cations, including necrotizing enterocolitis (NEC). NEC 
has a high mortality and this inflammatory condition 
affects up to 5-10% of the hospitalized preterm infants 
[6]. Nevertheless, the immature GIT is able to respond 
to the first enteral feeding and bacterial colonization with 
a large increase in tissue mass and digestive capacity 
[5,7,8]. However, the responses are highly diet-dependent 
and even a short period of feeding with suboptimal milk 
formulas may trigger the events that lead to NEC: mucosal 
hemorrhage, bacterial overgrowth, dysfuctional immune 
response, and failure of multiple organs [9]. Studies of the 
intestinal proteome suggest that the intestinal mal- 
adaptation in preterm neonates can be explained by a 
feeding-induced decrease in intestinal metabolism and 
stress response [8,10], but the molecular mechanisms 
remain unknown. 

Given the high rate of preterm birth (10-15%) world- 
wide [11], and the difficulties in feeding preterm infants 
[12], it is important to know the molecular mechanisms 
that guide pre- and postnatal GIT maturation in preterm 
neonates. Previous proteomic and transcriptomic studies 
demonstrated that preterm and term pigs differ markedly 
in the intestinal response to the first feeding [8] and both 
diet and bacterial colonization play a role [10,13]. Studies in 
pig fetuses and germ-free preterm pigs [5,7,13] showed that 
feeding-related responses occurred, even in the absence of 
bacterial colonization [14,15]. Together, these data suggest 
that there are genome-wide, postnatal feeding-dependent 
changes in the gene transcription of premature gut. Due to 
the essential role of epigenetics in developmental regulation 
of gene function, a genome-wide characterization of epi- 
genetic modifications in the preterm intestine will help to 
understand the mechanisms whereby the intestine of pre- 
term neonates adapts to its new environment. 

We hypothesized that epigenetic modification is a cen- 
tral mechanism whereby transcriptional changes occur 
in the intestine during the perinatal period of preterm 
newborns. Using our well-established newborn piglet 
model of intestinal adaptation and NEC, we performed a 
Reduced Representation Bisulfite Sequencing (RRBS) to as- 
sess the intestinal DNA methylation in four groups of pig- 
lets: 1) Full-term newborn (Od-term), 2) preterm newborn 
(Od-preterm), 3) healthy preterm newborns fed with infant 
formula (4d-preterm), and 4) preterm newborns fed with 
formula and showing NEC symptoms (4d-preterm-NEC). 
By group comparisons, we investigated the effects of 



advancing prenatal (PN) age (Od-term vs. Od-preterm), 
increasing neonatal (NN) age (Od-preterm vs. 4d- 
preterm) and NEC (4d-preterm vs. 4d-preterm-NEC) 
on DNA methylation of intestinal genes. This initial 
study showed substantial descendent change of DNA 
methylation in a genome-wide scale that was associ- 
ated with both PN and NN effects. Four genes 
(CYP2W1, GPR146, TOP1MT, CEND1) related to intes- 
tinal nutrient metabolism were significantly hyper- 
methylated in their promoter CpG islands (CGIs) and 
showed corresponding down-regulation of transcrip- 
tion levels, in healthy preterm newborns fed with in- 
fant formula for four days (4d-preterm). The observed 
promoter hypermethylation suggests profound epigen- 
etic effects of the first days after birth on the GIT in 
preterm pigs. Further studies are required to evaluate 
the observed methylation changes and how they vary 
with different diets and GIT bacterial colonization. 
The results may help to identify new tissue markers of 
GIT maturation and maladaptation-associated compli- 
cations, such as feeding-induced NEC. 

Results 

RRBS data generation and characteristics of porcine 
intestinal methylome 

To characterize the intestinal methylome of preterm 
piglets associated with the adaptation around birth and 
development of NEC, we collected empty, 2 cm full-wall 
sections from the mid intestine (50% intestinal length 
from the stomach) from four groups of newborn piglets, 
all delivered by caesarean section. Two groups were 
newborn, unfed pigs, either full-term (Od-term) or pre- 
term (Od-preterm, each n = 2) that were sacrificed within 
two hours after birth. The other two groups consisted of 
preterm piglets fed with infant formula for four days 
following the rearing and feeding protocol specified pre- 
viously [16], but showing severe symptoms of NEC in 
regions of the intestine other than the mid intestine that 
was selected for the analyses (4d-preterm-NEC) or being 
healthy (4d-preterm, each n = 3; Additional file 1: Table 
SI and S2). 

Genomic DNA of mid-intestines from each individual 
from the same group of piglets was extracted and pooled 
together to construct a single RRBS library, and in 
total four libraries were then sequenced using Illumina 
HiSeq2000 sequencer. RRBS was developed and widely 
used to measure DNA methylation of high-CG regions 
at single base-pair resolution [17,18]. In two previous 
studies, we thoroughly evaluated the genomic coverage 
as well as repeatability of this technology in human 
samples [19,20]. A dataset size of 5 Gb will ensure a 
good genomic coverage and coverage depth [19]. Fur- 
thermore, by correlating the methylation levels of 
commonly covered CpG sites from two independent 
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experiments, we gained a Pearson's coefficients of 0.98, 
which suggested for good technical repeatability of 
RRBS [20]. As this is the first time to use RRBS in the 
pig genome, we first performed in silico simulation on 
the genomic fragments that were enriched by digesting 
the pig genomic DNA with Mspl, selecting 40-220 bp 
fragments and performing a 50-bp pair-end sequen- 
cing. As a result, a maximum of 3.41 million distinct 
CpG dinucleotides (12.17% of all CpGs) were recov- 
ered from the sequenced fragments. Similar to other 
mammalian genomes [20], around half (52.30%) of the 
CGI CpGs in the porcine genome can be enriched by 
RRBS (Additional file 1: Table S3). We generated a 
total of 23.15 gigabases (Gb) clean sequencing reads 
from four libraries of pooled intestinal samples (5Gb 
minimal for each library), reaching a minimal coverage 
depth of 12.5 per strand for each (Additional file 1: 
Table S2). By mapping the clean reads to the reference 
genome, we obtained a similar distribution of sequen- 
cing reads/fragments in the four libraries, indicating a 
highly efficient enrichment of CpG dinuceotides and 
relatively unbiased RRBS library construction from our 
experimental procedure (Additional file 1: Table S3). 

After bisulfite treatment and PCR amplification, a meth- 
ylated cytosine is read as "C" while a non-methylated cyto- 
sine is read as "T" in the sequencing. We calculated the 
ratio of total "C" reads to total sequencing reads to initially 
determine global level of methylcytosines (mCs) for each li- 
brary with respect to the sequence context. We found that 
non-CG methylation was rare (<1%), compared with CG 
methylation (~53%), which is consistent with previous 
observations in mammalian somatic genomes (Additional 
file 1: Table S2). Besides, RRBS technology enriches high- 
CpG regions in the genome. Thus, we analyzed CpG 
methylation in the current study. Absolute methylation 
level for each cytosine in the library was defined as the ratio 
of "C" reads to total sequencing reads for a specific cyto- 
sine. We then investigated the DNA methylation level ac- 
cording to the chromosome features, and observed that the 
average methylation level correlated negatively with 
chromosome length (Pearson's r = -0.6, P = 7.8e - 05), re- 
peat density (r = -0.6, P = 7.4e - 19), GC content (r = -0.5, 
P = 4.7e-19) and CpG density (r=-0.2, P = 5.6e-20). 
There was no obvious correlation between methylation 
level and gene density (r =0.1, P =2.4e - 24) or observed 
versus expected number of CpG (CpGo/e) ratio (r=-0.0, 
P =4.6e - 19) (Additional file 2: Figure SI). The negative 
correlation between the methylation level of individual 
CpGs and GC content might be explained by the observa- 
tion that CpGs in the regions of high CpG density (so- 
called CpG islands, CGIs) [21] tended to be unmethylated, 
whereas CpGs in low-density regions tended to be methyl- 
ated, which is consistent with previous reports in mam- 
malian genomes [2]. These results indicated that porcine 



intestinal methylomes can be detected by RRBS technol- 
ogy with good representation, thereby ensuring accurate 
examination of DNA methylation alterations. 

Differential methylation patterns among four groups 
of piglets 

To examine the global pattern of the methylome, relative 
proportion of CpG sites, whose DNA methylation levels 
were classified into quintiles, were firstly depicted for 
each group of samples (Figure 1A). A different pattern 
of 5-mC distribution was observed between the preterm 
groups (4d-preterm, 4 d-preterm-NEC) and the newborn 
groups (Od-preterm, Od-term). Especially, the proportion 
of highly methylated cytosines (80-100%) decreased in 
the preterm groups. Hierarchical clustering analysis by 
"pvclust" algorithm [22] on methylation levels of all gen- 
omic CpG sites also revealed a clear separation between 
the 4 d neonatal and 0 d newborn piglets (Figure IB). 
Furthermore, the height of the "pvclust" tree diagram 
showed the intra-divergence within a cluster. Therefore, 
a higher height of the cluster with unfed, newborn 
groups (Od preterm and term) suggested more clear dif- 
ferences between these two groups than the postnatal 
groups (4 d preterm groups). Correspondingly, a propor- 
tional increase of low-methylated cytosines (0-20%) were 
observed between the Od-term and Od-preterm piglets, 
while the patterns of methylated cytosines between 4d- 
preterm and 4d-preterm-NEC pigs were more similar 
(Figure 1A). In addition, we observed that the unfed, 
newborn pig groups (Od) had markedly higher global 
methylation level than the two neonatal groups (4d) in 
gene bodies and downstreams. The methylation differ- 
ence was 3.4% and 1.8% in gene body and downstream 
5 k, respectively, while the difference in upstream 5 k 
was only 0.8% (Figure 1C). Considering the uneven 
distribution of methylation differences in different gene 
regions, and the similarly high bisulfite conversion rate 
among these libraries, as reflected by non-CpG methyla- 
tion levels (<1%), this result is not likely to be due to 
experimental bias. 

Next, we carried out pair-wise comparison to identify 
the differentially methylated regions (DMRs) through 
the whole genome using a window sliding strategy [23] 
(see Methods). X-chromosomes were excluded from the 
pairwise comparison considering the DNA hypermethy- 
lation that occur during X-chromosome inactivation 
[24]. There was a clear PN effect (Od-preterm vs. Od- 
term) on DNA methylation levels, and both these groups 
differed from healthy preterm piglets fed formula for 
4 days (4d-preterm). The DMRs between 4d-preterm 
and 4d-preterm-NEC pigs differed in their susceptibility 
to NEC lesions. As a result, the number of DMRs was 
3801 (averagely 117 bp in length, 14 CpGs) for Od- 
preterm vs. Od-term (PN effect), 5778 (averagely 122 bp 



Gao et al. BMC Genomics 2014, 15:716 
http://www.biomedcentral.com/1471-2164/15/716 



Page 4 of 14 



A 



B 



1.0 n 




0.8-1.0 
0.6-0.8 
0.4-0.6 
0.2-0.4 
0-0.2 



o 
o 



CD 
00 

— o 



CD 
X 



CM 
CO 
O 



00 
CN 
O 



Cluster dendrogram with AU/BP values (%) 
au bp 



100 100 



edge # 



100 100 



S 
E 

Q. 



1 



CD 

£ 

Q. 

-6 



o 

LU 

z: 



CD 



CD 
O 



CD 

E 
O 

Q. 
O 



o 
o 



Overview of methylation level in the genie regions. 



Distance: correlation 
Cluster method: average 



Od-term 
Od-preterm 
4d-preterm 
4d-preterm-NEC 




TSS up 5K 



Gene body 



2 3 4 
TTS down 5K 



Figure 1 Global DNA methylation profiles. (A) Global pattern of CpG methylation in the four groups. The methylation levels of all CpGs were 
calculated and categorized into five color-coded methylation states. CpG coverage (y-axis) shows the proportion of CpGs covered with different 
methylation levels. (B) The "pvclust" tree diagram clustering based on methylation of CpGs in the whole genome. (C) Average CpG methylation 
levels of different genetic regions. Five kilo base regions upstream and downstream of each gene were divided into 100-base pair (bp) intervals. 
Each coding sequence was divided into 20 intervals (5% per interval). 



in length, 13 CpGs) for Od-preterm versus 4d-preterm 
(NN effect), and 1259 (averagely 79 bp in length, 11 
CpGs) for 4d-preterm vs. 4d-preterm-NEC (NEC effect, 
Additional file 1: Table S4). More than 66% of the DMRs 
were distributed in intron and intergenic regions, while 



7.3-8.5% of the DMRs were located within upstream 2 k 
to TSS regions, indicating broad but similar distribution 
patterns across the whole genome (Figure 2A). We then 
cross-matched PN- and NN-DMRs to identify the DMRs 
that overlapped with at least half of the CpG sites in 
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each DMR. Interestingly, a Pearson correlation coefficient 
of 0.53 was obtained, showing a good correlation for most 
of the common CpG sites between PN- and NN-DMRs 
(Additional file 2: Figure S2). 

Decreased DNA methylation with different prenatal and 
neonatal age for genes affecting gut maturation 

The GIT normally undergoes a tremendous maturation 
in the perinatal period in order to cope with changes in 
the environment [5]. In order to investigate the intestinal 
DNA methylation alterations prior to normal birth, we 
compared the groups of unfed, newborn preterm and term 
pigs where the majority (76.9%) of the DMRs showed 
decreased methylation level in term piglets, compared with 



preterm piglets (Additional file 1: Table S4). Further 
changes occurred in the first few days after birth, as 
shown by the comparison of 4d vs. Od preterm pigs, 
where the majority of DMRs (79.6%) showed lowered 
methylation level in the 4d group. Among the CpG 
sites shared between the PN- and NN-DMRs, the ma- 
jority of the shared CpG sites showed concordance of 
decreasing methylation due to increasing gestational 
age at birth and the effect of the first few days after 
birth (Additional file 2: Figure S2). Thus, both prenatal 
age at birth and the first few days after birth have 
strong impacts on the methylome. 

As promoter hypermethylation is typically associated 
with the repression of gene transcription, we then focused 
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Figure 2 Summary of DMRs associated with prenatal maturation (PN-DMRs, Od-term vs. Od-preterm), neonatal response (NN-DMRs, 
4d-preterm vs. Od-preterm) and NEC development (NEC-DMRs, 4d-preterm-NEC vs. 4d-preterm). (A) The distribution of all DMRs identified 
with respect to prenatal maturation (PN-DMRs, 3801), neonatal response (NN-DMRs, 5778) and NEC (NEC-DMRs, 1259) in genetic regions. (B) Violin box 
plots show the average methylation levels of DMRs with respect to PN-DMRs, NN-DMRs and NEC-DMRs. The white points depict median values, the black 
boxes connect the 25th and 75th percentiles and the thin black lines connect the lower adjacent value to the upper adjacent value. The yellow area 
depicts a density trace, plotted symmetrically above and below the horizontal box plot. Student's t-test was used to identify significant differences of 
methylation levels of DMRs located in promoters between each pairs. (Q Venn diagrams show the result of the cross-matching genes with DMRs 
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on the genes that contained DMRs in their promoters. In 
total, 217 and 284 DMRs were located in the promoters of 
203 and 261 genes for the PN and NN comparisons 
(Additional file 1: Table S4). The majority of the DMRs 
in promoters (74.7% and 71.8%, respectively) showed 
decreased methylation levels in response to advancing 
prenatal and neonatal age (PN and NN effects, Figure 2B). 
A high proportion of PN-DMRs (101 out of 203) over- 
lapped with CGIs, corresponding to 91 genes, which were 
then subjected to KEGG (Kyoto Encyclopedia of Genes and 
Genomes) pathway and disease enrichment analysis using 
WebGeStalt [25,26]. Five KEGG pathways with an adjusted 
P value less than 0.05 were considered as significantly 
enriched. These included the following pathways: gly- 
cosaminoglycan biosynthesis, oxidative phosphoryl- 
ation, amyotrophic lateral sclerosis, complement and 
coagulation cascades and neuroactive ligand-receptor inter- 
action. Disease enrichment analysis indicated that malfunc- 
tion of these genes may lead to a broad spectrum of 
putative disease terms, such as pregnancy complications, 
hyperlipoproteinemia and type 2 diabetes (Additional file 1: 
Table S5). These results suggest that genes involved in in- 
testinal metabolism, innate immunity and environmental 
information processing undergo very marked DNA methy- 
lation alterations during the perinatal period, reflecting an 
important role of DNA methylation in regulation of GIT 
development [27]. 

For NN-DMRs (Od-preterm vs. 4d-preterm), 123 out of 
261 overlapped with CGIs, corresponding to 109 genes 
(Additional file 1: Table S4). To reveal the genes with DNA 
methylation alterations that were likely induced specifically 
by the environmental changes (including enteral feeding 
and gut colonization) during the first four days of neonatal 
life, rather than simple age effects alone, we filtered out 32 
genes (Figure 2C) that were commonly identified by cross- 
matching PN-DMR genes with NN-DMR genes and sub- 
jected the rest, 77 genes, to KEGG pathway and disease 
enrichment analysis. The KEGG pathway, focal adhesion, 
and adhesion and subarachnoid hemorrhage diseases, were 
enriched from these genes (Additional file 1: Table S5). 

In contrast to the DMRs associated with prenatal (PN) 
and neonatal (NN) age (see above), only around half 
(561) of the 1259 identified intestinal DMRs for NEC 
showed decreased methylation level in the 4d-preterm- 
NEC vs. the 4d-preterm piglets (Additional file 1: Table S4). 
In total, 49 NEC-DMRs were found located in pro- 
moters of 47 genes, 53% of which showed decreased 
methylation levels in the NEC piglets. Thus, about equal 
proportion of intestinal methylation changes, increasing 
or decreasing, were observed in response to NEC (Figure 2B 
and 2C). The 25 genes containing 27 NEC-DMRs that 
overlapped with CGIs. Significantly, four KEGG pathways 
were enriched out of these 25 genes (P value < 0.05), includ- 
ing calcium signaling pathway, neuroactive ligand-receptor 



interaction, oxidative phosphorylation and alzheimer's dis- 
ease. Disease analysis showed an enrichment associated 
with infection and immune disease, like asthma (Additional 
file 1: Table S5). Interestingly, ten out of the 25 genes 
overlapped either with PN-DMR or with NN-DMR 
genes (Figure 2C), suggesting methylation of these 
genes was highly related to NEC outbreak, advancing 
prenatal age and advancing neonatal age. 

Development of BSP combined with high-throughput 
sequencing for methylation validation 

As we only used 2-3 piglets in each group for initial 
genome-wide screening of methylation alterations, fur- 
ther validation of DNA methylation alteration of identi- 
fied key genes were performed on larger sample sets. To 
examine the methylation of one single gene, bisulfite- 
sequencing PCR (BSP) combined with molecular cloning 
of the amplified fragments and conventional Sanger se- 
quencing is the most routinely used method. However, 
the cost of deep sequencing of a large set of genes is 
high. Furthermore, greater variability in the methylation 
of CpG sites has been observed than in the direct se- 
quencing strategy, like pyrosequencing [28]. Therefore, 
we developed a new BSP-based procedure by using 
high-throughput Illumina HiSeq sequencing (HiSeq- 
BSP) instead of Sanger cloning sequencing (Sanger- 
BSP). To do so, PCR products from multiple genes of 
one sample were pooled together for library construc- 
tion by Illumina Pair-End protocol, and different li- 
braries were barcoded for high-throughput sequencing 
using HiSeq2000 (see Methods). A deep sequencing of 
the PCR products can thereby be achieved for each 
PCR fragment cost-effectively. 

To evaluate the performance of this new method, we 
first applied both HiSeq-BSP and Sanger-BSP strategies 
on seven genes. Based on the uniquely mapped HiSeq 
sequencing reads, 80% of CpG sites on the PCR se- 
quences could be covered with an average of 1000 fold 
sequencing depth (Additional file 1: Table S6), while an 
average of 20 clones were generated for each of the PCR 
fragments in Sanger-BSP. A scatter plot depicting methy- 
lation levels of all CpG sites showed a general consistency 
of methylation levels of CpG sites in 5 out of 7 genes be- 
tween these two approaches (Figure 3A and 3B), indicat- 
ing a high feasibility of the new HiSeq-BSP approach. 
However, Sanger-BSP tended to produce higher methyla- 
tion levels than HiSeq-BSP, especially for genes, CDX1 
and MEST (Figure 3A and 3C). We further performed a 
Sanger-BSP on CDX1 in the same samples independently. 
We found that a second cloning procedure generated, not 
only a higher, but also a different methylation profile for 
the same gene (Figure 3C). Twenty clones might not be 
deep enough to secure a random selection of the PCR 
fragments, and the additional cloning of bisulfite-treated 
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DNA into bacterial expression vectors might increase its 
variability of methylation detection. These results indicate 
that HiSeq-BSP had better accuracy than traditional 
Sanger-BSP, due to high-depth direct sequencing. 

Expression of key genes repressed by hypermethylation 
of CGI 

We then collected mid-intestinal samples from another 24 
piglets (each group n = 6) for further validation (Additional 
file 1: Table SI). In order to address the issue of potential 
inter-individual variations within the same group, we chose 
to examine the methylation level of key genes in each pig- 
let, instead of a pooling strategy used in high-throughput 
RRBS analysis. We selected DMRs for further validations 
based on two criteria: First, with high co-methylation varia- 
tions (difference value of methylation levels >20%), second 
with at least 5 CpGs and length > 30 bp. As a result, 19 
PN-DMRs (18 genes), 15 NN-DMRs (15 genes) and 8 
NEC-DMRs (6 genes) were successfully amplified from 
each of the 24 piglets and subjected to HiSeq-BSP sequen- 
cing. Considering that inter-individual variations might lead 
to a shift of DMR regions, we re-performed pair-wise 
comparison to search for DMRs on the amplified PCR 
fragments between two groups, 6 piglets for each group 
(Additional file 1: Table S6). Four (22%) PN-DMR genes, 
9 (60%) NN-DMR genes and 2 (33%) NEC-DMR genes 
were consistently discovered with DMRs in the validation 
results, which all agreed with RRBS results (Additional file 1: 
Table S6). We calculated coefficients of variation (CV) for 
each CpG site among six piglets from each group, in order 
to determine the inter-individual variations. The CpG sites 
with high CVs were located within ten genes, nine of 
which failed the validation (Additional file 1: Table S7). 
The highest percentage of validation was achieved in 
NN-DMR genes, which might suggest that the most stable 
pattern of methylation changes is seen in response to the 
first 4 days of neonatal life in preterm neonates. 

We then tested whether the methylation changes affected 
the transcription level of the subset of genes with validated 
DMRs. Eleven genes were successfully subjected to gene 
expression test using quantitative RT-PCR, including three 
PN-DMR genes {DEGS2, ADAM30 and KLHL35), six 
NN-DMR genes {DLEU7, CEND1, RAB11FIP3, CYP2W1, 
GPR146, TOP1MT) and two NEC-DMR genes (CDH22 
and ZPBP). These genes are implicated in a variety of 
biological processes related to organ development and 
intestinal functions. As a result, four NN-DMR genes 
(CYP2W1, GPR146, TOP1MT and CEND1) were sig- 
nificantly down-regulated at the transcription level 
(Od-preterm vs. 4d-preterm, fold change > 2, P value < 
0.01, Figure 4), which were negatively correlated with 
their hypermethylation in promoter CGIs. CYP2W1 
[29], GPR146 [30] and TOP1MT [31] are important for 
nutrient metabolism in the intestine. CEND1 plays a 



role in neuronal differentiation [32], thus in development of 
the enteric nervous system. The down-regulation of these 
genes due to promoter hypermethylation may significantly 
affect intestinal metabolism and development, which may 
help to explain the difficulties of preterm neonates to toler- 
ate enteral feeding and the down-regulation of proteins in- 
volved in intestinal metabolism and motility [8]. Among 
the three tested PN-DMR genes and two NEC-DMR genes, 
no significant difference was observed in their expression 
levels between the groups (Figure 4), suggesting that 
the effect of increasing gestational age (PN effect) and 
development of NEC on the expression of these genes 
was limited. 

Discussion 

GIT maturation involves a continuous cascade of growth, 
differentiation, and renewal of epithelial cells [33]. 
Rapid GIT adaptive changes occur at the phenotypic 
and genome-wide transcription levels during the first 
few days of life in both preterm and term neonates [5]. 
It is only in preterm neonates, however, that the initial 
exposure to neonatal life may lead to intestinal dysfunc- 
tion and NEC [5]. Previous studies in both preterm pigs 
and infants indicate that the intestinal phenotype sensi- 
tive to develop NEC is characterized by maladaptation 
of the intestine in response to postnatal conditions, as 
shown by impaired intestinal metabolism, perfusion, 
motility and bacterial defenses [5]. Environmental fac- 
tors, coupled with impaired host-response mechanisms, 
appear to determine NEC sensitivity, and neither epidemio- 
logical studies in infants, nor the numerous studies in 
preterm pigs, suggest a large component determined by 
maternal or paternal genetics. Interestingly, it is mainly 
when infant formula, not mother's milk, is used as the first 
enteral diet that the intestine fails to adapt. It is not known 
how the perinatal adaptation of the immature intestine in- 
volves changes at the gene level, and to which extent these 
changes occur prenatally or postnatally along with the 
introduction of milk feeding and bacterial colonization. We 
now show that the immature intestine undergoes rapid 
changes in gene methylation during the last 10 d of gesta- 
tion, but even more dramatic changes occur when the im- 
mature intestine is exposed to just 4 d of formula-feeding 
during postnatal life. Neonatal hypermethylation of genes 
related to intestinal metabolism and enteric nervous system 
may partly contribute to the inappropriate response of im- 
mature GIT to feeding and bacterial colonization. The 
hypermethylation may potentially cause long-term prob- 
lems of intestinal functions in preterm neonates although 
the effects could also be more temporary, as DNA methyla- 
tion is highly dynamic throughout the lifespan of mammals 
[3,34,35] . Further studies on the long term effects of differ- 
ent neonatal feeding patterns may provide an answer to this 
important question. 
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Figure 4 Neonatal promoter hypermethylation associated with gene transcription repression. After HiSeq-BSP, genes with consistent 
DMRs with RRBS results were subjected to real-time qPCR analysis to determine the differences of their expression levels with respect to prenatal 
age (PN), neonatal response (NN) and NEC disease (NEC). Four NN-genes (CYP2W1, GPR146, TOP1MT, CEND1), related to intestinal metabolism, 
were significantly hypermethylated in their promoter CGIs (left), and transcriptionally downregulated in the 4d-preterm compared with Od-preterm (right). 
The quantitative ratios are normalized to the expression of GAPDH. The data are representative of three similar experiments and are displayed as the 
mean ± SD. * P < 0.01 as evaluated using the Student's t-test and with a fold change > =2. 



We studied for the first time DNA methylation changes 
in a mammalian immature intestine around the time of 
birth. Whole tissue analyses without separating different 
tissue layers or cell types were performed in this study 



because developmental programs are often a combined re- 
sult of interactions among many tissue layers and cell types. 
Thus, studying only one single cell type would not provide 
information about the changes at the whole tissue level, 
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although the variability in results would have been less if 
just one cell type had been studied. Similar studies in pre- 
term infants are impossible due to the inaccessibility of 
tissue samples and piglets are probably the best model to 
study early GIT development and disease [36] and thereby 
identify preventive strategies that help to avoid intestinal 
disease in preterm infants [37]. The results of our study 
highlight the importance of intestinal gene methylation 
during the last 10 days of gestation and the first 4 days after 
preterm birth, and identified several key intestinal genes 
that were affected. During both of these critical adapta- 
tion periods, intestinal DNA methylation increased in 
a genome-wide scale, indicating the repression of cor- 
responding genes, and that this age-related methyla- 
tion program was at least partly independent of the 
environmental changes at birth. On the other hand, the 
marked effects of the first 4 days of life (NN-DMRs 
relative to PN-DMRs) suggest strong environmental 
effects in addition to the predefined ontogenetic pro- 
gram. From our results, we cannot know whether the 
key neonatal adaptation factor is formula feeding, bac- 
terial colonization or some other postnatal factors, but 
some previous results suggest that enteral feeding may 
be central for the changes. Enteral feeding alone, with- 
out bacterial colonization, can induce enormous meta- 
bolic changes in the immature GIT as germ-free fetal 
or preterm newborn pigs also show marked maturation 
after receiving enteral feeding for just 24 h [7,38]. 
Intestinal dysfunction, leading to NEC, depends on 
bacterial colonization, but in this study we observed 
limited gene methylation changes in response to devel- 
opment of NEC (NEC-DMRs). For analyses in severely 
NEC-affected pigs, we deliberately selected individuals 
with a mid intestinal section that appeared healthy at 
tissue collection. This was to avoid that conclusions re- 
garding methylation would not merely reflect the 
process of tissue pathology and necrosis. The results 
show that high sensitivity to NEC did not in itself in- 
duce profound changes in tissue methylation and we 
suggest that formula feeding itself is the most deter- 
mining factor for intestinal gene methylation changes 
just after preterm birth. More studies on different feed- 
ing modes, diets and degrees of gut colonization are 
required to know whether intestinal gene methylation 
is affected by type of diet, just like NEC sensitivity. An 
optimal neonatal diet may promote long term health of 
preterm neonates by inducing permanent changes to 
the intestinal methylome. 

The advancing gestational age at birth resulted in PN- 
DMRs in promoters associated with genes enriched in 
KEGG pathways that are considered important for gut 
development. "Oxidative phosphorylation" is a meta- 
bolic pathway that is indicative of mitochondrial activ- 
ity and central for intestinal metabolism [39,40], while 



"glycosaminoglycan biosynthesis" is key for the intes- 
tinal structure and first line defense against pathogens, 
and they are down-regulated in intestines affected by 
NEC [41,42]. It appears that glycosaminoglycan bio- 
synthesis is regulated by methylation changes, and 
after birth, mucin synthesis was also down-regulated in 
formula-fed pigs compared with colostrum-fed pigs [43]. 
The "Complement and coagulation cascades" pathway is a 
part of innate immunity as a nonspecific defense against 
pathogens [44], suggesting epigenetic changes regulates 
this protective mechanism during the prenatal period. 

We observed a relatively low rate of DMR validation 
(4/15). The clear inverse relationship between DMR level 
and gene transcription (0/3) for a larger group size of 
samples from newborn preterm and term animals (PN 
effect) suggests high individual variation among piglets. 
High individual variation could also explain the relatively 
small average value of methylation differences in the key 
KEGG pathway genes between these two groups. Re- 
gardless, it is possible that the immaturity-associated 
DNA methylation pattern may lead to unsuccessful 
activation or silencing of key genes involved in normal 
prenatal gut development, thus helping to explain the 
higher risk for NEC in preterm neonates [5,45]. For ex- 
ample, we observed promoter hypermethylation, and 
lowered expression, of the DEGS2 gene in the newborn, 
preterm piglets. This gene plays an essential role in the 
sphingolipid synthesis in intestine [46], which in the neo- 
natal period may be supported by milk sphingolipids [47]. 

In response to the first few days of neonatal life, the 
majority of NN-DMRs showed hypomethylation (4d- 
preterm vs. Od-preterm), which is consistent with a gen- 
eral increase in the gene functions to support intestinal 
growth and functions after birth. Four hyper-methylated 
NN-DMR genes were identified in their promoter CGIs 
of the 4d-preterm pigs, and the transcription of these 
genes was correspondingly down-regulated. Three of 
these genes are related to intestinal metabolism. GPR146 
is a G protein-coupled receptor, which acts as nutrient 
chemosensor in the intestine [48] and transmits proinsu- 
lin C-peptide signaling [49]. CYP2W1, an enzyme of the 
cytochrome P450 superfamily, catalyzes the reactions in- 
volved in drug metabolism and synthesis of cholesterol, 
steroids and other lipids [29]. TOP1MT, a mitochondrial 
DNA topoisomerase, is critical for mitochondrial integrity 
and cellular energy metabolism [31]. Decreased formula- 
induced expression of the above metabolic genes may 
disturb normal intestinal nutrient metabolism and may 
thereby predispose to complications. In accordance 
with these observations, we previously found that a 
switch from parenteral to enteral feeding with formula 
(but not mother's milk) rapidly induces diet-dependent 
histopathological, functional, and proinflammatory changes 
in the immature intestine [9], as well as decreased 
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expression of specific proteins involved in intestinal metab- 
olism [8]. We now provide evidence that these immediate 
neonatal responses may be mediated by more fundamental 
changes to the genes involved, and that this may have long 
term consequences for intestinal functions in preterm neo- 
nates. Further studies on older preterm and term piglets, 
fed with different diets in the immediate postnatal period, 
can provide answers to these hypotheses. Potentially, the 
feeding time, dose, diet (infant formula or mother's milk) 
and microbial colonization could exert differential effects 
on intestinal gene methylation. Dietary methyl donors may 
also be important, as shown by the observed effect of folate 
on DNA methylation [4]. Infant formula is very different 
from mother's milk in the contents of bioactive nutrients 
and immunological, antimicrobial and growth-stimulating 
factors [50]. 

Development of NEC in preterm infants is characterized 
by a severe inflammatory response, including epithelial 
degeneration, edema and mucosal destruction. Our current 
study demonstrated that the DNA methylation changes as- 
sociated with NEC were enriched in the "focal adhesion" 
KEGG pathway and relevant to "adhesion" and "subarach- 
noid hemorrhage", suggesting that immature DNA methy- 
lation might be related to NEC. However, the pairwise 
comparison between healthy and NEC preterm pigs only 
revealed two genes with notable changes in DNA methy- 
lation levels and we could not document corresponding 
changes at the corresponding gene expression levels. 
CDH22 is a cadherin protein associated with colorectal 
cancer disease [51], while association of ZPBP with 
inflammatory bowel disease was shown in a previous 
genome-wide association study [52]. Despite the marked 
phenotypic changes in the intestine related to occurrence 
of NEC, no marked changes in intestinal DNA methyla- 
tion was found, at least not within the first days of neo- 
natal life. 

Conclusions 

In conclusion, our data revealed marked developmental 
changes in DNA methylation in the immature prenatal 
intestine and especially during the first days after preterm 
birth. The marked structural and functional changes that 
occur in the intestine at this critical time were accompanied 
by large changes in the intestinal methylome, which likely 
induce long-lasting effects beyond the neonatal period. We 
speculate that enteral feeding alone can bring profound 
changes to DNA methylation levels for intestinal genes. 
Diet- and feeding-related alterations in the DNA methyla- 
tion of specific genes may contribute to the high sensitivity 
to NEC. Further studies on intestinal DNA methylation 
changes, both short and long term, in response to neonatal 
feeding and bacterial colonization are needed. These could 
provide important information about the mechanisms of 
intestinal adaptation in all newborns, and potential targets 



for interventions against gut diseases in preterm newborns. 
Preterm birth is associated with long term complications in 
many organs and it is important to determine whether peri- 
natal methylation changes can be manipulated by optimiz- 
ing key environmental factors, such as diet and bacterial 
colonization at this critical time. 

Methods 

Animal treatment and tissue collection 

Six term pigs were delivered close to full term (115-116 
days gestation, by cesarean section) from two pregnant 
sows (Large White x Duroc, Askelygaard, Denmark), serv- 
ing as the Od-term group. Within 6 hours after birth, these 
pigs were euthanized (sodium pentobarbitone, 200 mg/kg, 
intravenous), and an approximately 2-cm-long intestinal 
segment left to the middle point of whole small intestine 
was taken from each pig, snap frozen in liquid nitrogen and 
stored at -80°C This method was adopted in the following 
samplings of intestine. Under anesthesia, 18 preterm pigs 
were obtained from five pregnant sows (Large White x 
Duroc; Askelygaard, Denmark) by cesarean delivery at 90- 
92% (104-107 days) gestation, as described earlier [38]. Six 
of these pigs were euthanized within 6 hours after birth 
without any feeding (Od-preterm group). The intestine seg- 
ments were collected and frozen at -80°C. The remaining 
12 preterm pigs were reared individually in infant incuba- 
tors (Air-Shields, Hatboro, PN) and fitted with orogastric 
catheters for parenteral nutrition and orogastric feeding 
tube. The pigs were given parenteral nutrition at a rate of 
4 mL/(kg x h) on dl and 6 mL/(kg x h) on d2 as well as 
minimal enteral nutrition: 24 mL/(kg x d) on dl and 
40 mL/(kg x d) on d2 with infant formula. The parenteral 
nutrition was based on a commercially available product 
(Kabiven, Fresenius Kabi) and adjusted in nutrient compos- 
ition to meet the requirement of pigs. On d3, they were 
transferred to full enteral nutrition [120 mL/(kg x d)] for 
two more days before sacrifice. The mid intestine was also 
dissected out and frozen at -80°C. Detailed feeding 
protocol was previously described [16]. Pigs were con- 
tinuously monitored for clinical symptoms of NEC and 
at tissue collection the severity of NEC-like intestinal 
lesions was evaluated using a scoring system previously 
described [16]. Six pigs remained healthy at the time of 
tissue collection on day 4 (4d-preterm) while six other 
pigs had significant NEC lesions in their intestines 
(4d-preterm-NEC). For validation experiments, mid 
intestine samples from other 6 piglets for both groups 
were collected, which were treated in the same way in 
RRBS study. Information on all the samples is included 
in Additional file 1: Table SI. All the procedures on 
animals were approved by the National Committee 
on Animal Experimentation in Denmark (permit no. 
2012-15-2934-00193). 
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The RRBS library construction 

For each intestine sample, total DNA was prepared by 
proteinase K/phenol extraction and examined by agarose 
gel electrophoresis. Then, the extracted DNA of all in- 
testine samples within each group were pooled together 
for RRBS library construction, thereby four pools of 
intestinal genomic DNA were prepared for the four groups 
of piglets (Additional file 1: Table SI), respectively, and used 
for RRBS library construction as previously described [19]. 
Briefly, four ug of genomic DNA was digested with 100U of 
Mspl enzymes (NEB) in 100 uL reactions at 37°C overnight. 
After purification, the digested products were progressed 
by blunt-ending, dATP addition and methylated-adapter 
ligation. Two ranges of 160-240 bp and 240-340 bp frac- 
tions were excised from a 2% agarose gel, respectively. Then 
bisulfite conversion was conducted by ZYMO EZ DNA 
Methylation-Gold Kit™ (ZYMO) following the manufac- 
turer's instructions. The converted DNAs were amplified 
using JumpStart™ Taq DNA Polymerase (Sigma) by 11-13 
cycles (11 cycles for 160-240 bp and 13 cycles for 240- 
340 bp) of 94°C for 30 sec, 58°C for 30 sec, 72°C for 30 sec. 
RRBS libraries were then analyzed by Agilent 2100 Bioana- 
lyzer (Agilent Technologies) and quantified by PCR. 

RRBS sequencing and data processing 

The RRBS libraries were sequenced using Illumina 
HiSeq2000 Analyzer according to the manufacturer's 
instructions. Raw sequencing data was processed by 
the Illumina base-calling pipeline. Low-quality reads 
that contained more than 30% 'N's or over 10% of the 
sequence with low quality value (quality value <20) per 
read were omitted from the data analysis. The pig 
reference genome (Sscrofal0.2) by Swine Genome 
Sequencing Consortium (SGSC) was downloaded from 
Genbank database (ftp://ftp.ncbi.nih.gov/genbank/genomes/ 
Eukaryotes/vertebrates_mammals/ Sus_scrofa/ Sscrofal0.2/) . 
The annotation was downloaded from Ensemble (Sus_ 
scrofa.Sscrofal0.2.67.gtf.gz). The clean reads were aligned 
to the reference genome using SOAP aligner (Version 
2.01) [53] according to a previously published method 
[23] . The uniquely aligned reads that contained Mspl en- 
zyme digestion site in the ends were used for further 
analysis. 

Identification of DMRs, KEGG pathway and disease 
enrichment analysis 

To identify putative DMRs, a sliding window strategy 
was used as previously described with minor adjust- 
ments [23]: First, commonly covered CpG sites with se- 
quencing depth > =5X between two groups were selected 
as candidate sites. After bisulfite treatment and PCR 
amplification, a cytosine will be read either as "T" if it's 
unmethylated or "C" if it's methylated, following binomial 
distribution, as previously suggested [54]. Methylation level 



of individual cytosine can then be defined as the ratio of 
"C" counts to total counts of "C" and "T" in the sequenced 
reads for each individual cytosine. Therefore, a two-tailed 
Fisher's Exact Test was first used to test the "C" and "T" 
counts for each cytosine between two groups. Then, the 
first CpG with significantly differential methylation (auto- 
matically corrected P-value < 0.05) was used as an initial 
locus of DMR, and following candidate sites were merged 
into a candidate DMR according to following criteria: First, 
the distance between two neighboring candidate sites < = 
300 bp, in order to cope with Msp/-digestion based RRBS 
strategy. Second, all candidate CpG sites in the candidate 
DMR maintain the same methylation direction (hyper- or 
hypo-). Third, a candidate DMR must harbor 5 or more 
candidate CpG sites, among which at least 3 CpG sites 
should have Fisher's exact test P-value less than 0.05. Last, 
for each of above candidate DMRs, a Fisher's exact test was 
performed again based on the mean "C" and "T" counts for 
all the CpG sites within the candidate DMRs. A false dis- 
covery rate (FDR) adjustment was then performed on the 
P-values, using the R package of "P.adjust" that is based on 
BH method [55]. DMRs with Fisher test FDR P-value <0.05 
and difference of mean methylation levels between two 
sample < 20% were taken as statistically significant. KEGG 
pathway and disease enrichment analysis was performed 
using the WebGestalt web server (http://bioinfo.vanderbilt. 
edu/webgestalt/) [25,26]. 

Validation by bisulfite sequencing PCR combined with 
cloning and Sanger sequencing 

PCR primers were designed using the online MethPri- 
mer software (www.urogene.org/methprimer/index.html) 
and listed in Additional file 1: Table S8. Genomic 
DNA (500 ng) was converted using ZYMO EZ DNA 
Methylation-Gold Kit™ (ZYMO) and one-tenth of the 
elution products were used as templates. PCR amplifi- 
cation was carried out with a thermal cycling program 
of 94°C for 1 min, 30 cycles of 94°C for 20 sec, 50 ~ 60° 
C for 30 sec, 72°C for 40 sec and then final 4 min incu- 
bation at 72°C. For Sanger sequencing, PCR products 
were purified using the QIAquick Gel Extraction Kit 
(Qiagen) and subcloned, 20 colonies from each prod- 
uct were sequenced using the 3730 Genetic Analyzer 
(Applied Biosystems). 

Validation by bisulfite sequencing PCR combined with 
HiSeq sequencing 

PCR primers were designed using the online MethPri- 
mer software (www.urogene.org/methprimer/index.html) 
and listed in Additional file 1: Table S8. For Hi-Seq se- 
quencing, PCR products of multiple genes from one 
sample were quantified by Tangen gel system and pooled 
together equally. One ug of pooled PCR products was 
fragmented by a Covarias sonication system to a mean 
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size of approximately 200 bp followed by end-polishing, 
A-tailing, adapter ligation and 5-cycle PCR to generate a 
barcoded library. Barcoded libraries from different samples 
were then pooled together equally and used for cluster gen- 
eration and standard pair-end sequencing with 50 bp reads 
(PE50) using Illumina HiSeq 2000. After sequencing, the 
Illumina reads were post-processed and aligned to the pig 
reference regions (all PCR regions) using SOAP aligner 
(Version 2.01) [53] according to a previously published 
method [23] with default parameters that excluded reads 
with more than five mismatched bases. Multiple reads 
mapping to the same position were counted only once to 
remove potential bias from PCR. 

Quantitative real time PCR 

Total RNA was isolated using Qiazol and RNeasy Midi 
kit (Qiagen), and cDNA was synthesized using 500 ng 
total RNA in a total of 20 ul reaction mixture using the 
QuantiTect Reverse Transcription kit (Qiagen) accord- 
ing to the manufacturer's instructions. Primers for real- 
time quantitative PCR (RT-qPCR) were designed using 
Primer 5 software and listed in Additional file 1: Table S8. 
Real-time qPCR analysis was performed on ABI Prism 
7700 (Applied Biosystems, Tokyo, Japan) using SYBR 
Green real-time PCR master mix (Toyobo Co., Japan). 
Relative expression levels of objective mRNAs were 
calculated using the AACt method and normalized to 
GAPDH. All data are presented as the mean values ± 
S.E. Comparisons were made using the Student's t-test 
and a two-sided p-value < 0.05 was considered to indi- 
cate statistical significance. 

Data availability 

All RRBS sequencing and processed data was deposited 
in the Gene Expression Omnibus (GEO) with accession 
GSE53747. Validation data was deposited in GigaDB, 
the GigaScience database, with the unique identifier 
doi:10.5524/100043 [56]. 

Additional files 

Additional file 1: This file contains Table S1-S8. 
Additional file 2: This file contains Figure S1-S2. 
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